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ri '• 1 Introduction 

In two recent papers , Q] , we have developed an algorithm for the solution of linear differential systems 
of the type 

> 

)0. Z'(x)=p(x){D + R(x)}Z(x) (1.1) 

O , 

on an interval [X, oo). Here Z is an n-component vector, p is a scalar factor, D is a constant diagonal 

OC • 

matrix and R is a perturbation matrix such that R(x) = 0(x ) as x — + cxd, for some <5 > 0. The algorithm 



whose perturbation matrices are of successively smaller orders of magnitude as x — > oo. When the 

> 

X 

of the final system in the process, the solution also possessing that accuracy. The corresponding solution 



implements a repeated transformation process by means of which ( [l. l[ ) is transformed into other systems 



perturbation reaches a prescribed accuracy, the Levinson asymptotic theorem p] , H provides the solution 



of (1. 1) is then obtained by transforming back. 

The algorithm has two main aspects. The first is the algebraic one of generating symbolically the 
matrix terms which appear in the transformation process. In ^| and the algorithm was set up in such 
a way that this aspect is implemented in the symbolic algebra system Mathematica. The other aspect 
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is the numerical one of deriving from the algebra the values of the solutions of ( |l. l| ), given p, D and R. 
Included here are the determination of explicit error-bounds and procedures for handling large numbers 
of matrix terms. 

In this paper we develop these ideas in two different but not unrelated directions. First, we replace 
the constant matrix D by one which has entries of differing orders of magnitude as x — > oo. Our methods 
are sufficiently indicated by two orders of magnitude, and thus we consider the system 

Z'(x) = p(x){A(x) + R(x)}Z(x) (1. 2) 



where 

^ XD ^ 



A 



(1. 3) 



with constant and diagonal D and D, while A is a scalar factor such that X(x) — > oo as x — > oo. The 
amplification of our previous algorithm to cover this new situation is dealt with in sections 2-4 below. 
At the end of section 3, we state the Levinson asymptotic theorem in the form that we need, with the 
necessary accuracy incorporated. As in and a typical p(x) is x a (a > —1). However, the value 
a = — 1 also occurs in applications, and our second extension is to include a = — 1 in the scope of the 
algorithm. This matter is the subject of section 5. 

Both these extensions occur together in a system (1.2) which is associated with the generalised hy- 
pergeometric equation. Since this system can be solved in terms of special functions, it provides an 
independent check on the effectiveness of our algorithm, and this matter is discussed in sections 6 and 
7. 
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2 Large scalar factor A 

The system Ql. | ) is the first in a repeated transformation process and therefore we write it as 

Z[(x) = p(x){A 1 (x) + Ri(x)}Zi(x). 



In order to increase the potential for applications, we broaden (1. 3) to 

Ai(i) = A (x) +$i(x), 

where 

An = 



XD ^ 


, $1 = 


( Ai ] 


V D ) 




^0 A, J 



Here 

D = dg(di, ...,d N ), D = dg(d N+ i, ...,d n ) 
are constant matrices, each one having distinct entries, and 

X(x) — > oo (x — > oo). 

Further Ai and Ai are diagonal with 

Ai=o(A), A x =o(l) (s^oo). 



(2. 1) 



(2. 2) 



(2. 3) 



(2. 4) 



The situation considered previously in |^] and (lj is where N — and, in extending this earlier work, we 
shall be brief when there are similarities. 

The nature of R\ is as introduced in J|, section 4] and 0, section 3], that is, R\ is the sum of terms 
of differing orders of magnitude in the form 



Ri = Vu + V 21 



(2. 5) 



where 



and 



V n =o{l), V kl =o(V 3l ) (k>j) 



Ei=o{V^) 



(2. 6) 



(2. 7) 



as x — > oo. We arrange that 

dgV n = (2. 8) 

by transferring any non-zero diagonal terms in V\\ to $i. In (2. 6) and ( |2. 7| ), £i has a prescribed 
accuracy as x — > oo, while the Vji do not have that accuracy and they are successively replaced by 
smaller-order terms as we implement the transformation process which we discuss now. 
The sequence of transformations has the form 



{I + P m )Z m+1 (m = l,2,...) 



(2. 9) 



where P m = o(l) (x — > oo) and the definition of P m in (2. 15) below is an extension of the one given for 
N = in 1^, (4.11)]. Starting with ( |2. l| ), a typical transformed system arising from ( 2. 9] ) is 



Z m — p(A m + R rn )Z n 



(2. 10) 



where, corresponding to ( 2. 5)-( 2. 8) 



Rm — Vim + V2r 



(2. 11) 



with 



V km = 0(Vj m ) (fc > j), E m = 0{V^ m ) 

as x — > oo, and dg V\ rn = 0. Here /i will depend on m, and the process ends at m = M when ( |2. ll[ ) 
reduces to Rm — Em with Rm thus achieving the prescribed accuracy. 



In ( |2. 10 ) we write 



A m = A + $ ro , $ m = dg(A m , A m ) 



(2. 12) 



where, as we shall see, A TO and A m are obtained respectively from Ai and Ai by adding to them terms 
of smaller orders of magnitude. Thus 



A ro = A 1 + o(A 1 ), A m = Ai + o(A!) 



(2. 13) 



as x — > oo. The transformation ( [2. 9|) takes ( 2. 10 ) into 

Z r ' n+1 = p[A rn + + P^i-p^P^ + Vl m + A,nP,n - PrnKn + {Rrn - Vl m ) + R m P m }]Z m+1 . (2. 14) 

We have to define P m so that the dominant matrix Vi m is removed from ( |2. 14j ) and, at the same time, 
P m has a form which can be readily included in an algorithm. 

Dropping certain subscripts for clarity, we define the entries py in P m in terms of the entries i>ij in 
Vim by 



Pa = S 



Vij/{X(dj - d,)} (1 < i < N, 1 < j < N) 

-Vij/(Xdi) (1 < i < N, N + 1 < j < n) 

Vij/(Xdj) (N + 1 < i < n, 1 < j < N) 

Vij/(dj -di) (N + 1 < i < n, N + 1 < j < n) 



where i ^ j, and dg P m — 0. Then, by (2. 2) and (2. 12), we have 



Vim A-rnPm Pm-^m U m -\- & m P nl P m <$} Y 



where the entries Uij in U m are 



\- 1 (d j /d i )v lj (l<i<N, N + 1 < j < n) 
\- 1 {d i /d j )v ij (N + l<i<n, 1 < j < N) 



(2. 15) 



(2. 16) 



(2. 17) 



and zero otherwise. We note that U m has a factor A l , making U m — o(V\ m ) as required. To assess the 



size of the terms involving <£> m in ( 2. 16), we first write 



(2. 18) 



where Q m contains the entries involving A in ( 2. 15 ) and Q m contains the remaining entries Vij/(dj — di). 
Then, by ( |2. 16[ ) and ( |2. 12j ), 



^ rn Pra Pm^m — P m H~ Pm ~\~ Pjni 



where 



P m — 



V 



( a ^ 

— * rn 



J 



T — 

± rn. — 



P \ 



and 



Qm Qm 

V 

\ / 

Qm 

V 



o 



A, 



\ 
) 



(2. 19) 



(2. 20) 



Qm Qr 



V 



(2. 21) 



/ 



Again, it follows from ( |2. 12| ), ( g 13| ), ( |2. 15| ) and ( p. | ) that T m ,f m and T m are all o{V\ m ). Then, 
provided also that 

p- l P' m = o(V lw ), (2. 22) 

we can proceed as in |, (4.7)] and @, (2.13)] to express (7 + P m ) 1 as a geometric series in ( [2. 14| ) and 
then collate terms of the same order of magnitude to obtain ( 2. 10 ) - ( [2. ll[ ) with m + 1 in place of 
If the dominant term is initially denoted by SVn+i, we arrange that dg Vi im +i = by defining 



rn . 



A m + i — A m + dgS* m +i, Vl m+l — S'rjj + l — dgSVn+1- 



In particular, A m is built up from Ai by adding terms of successively smaller orders of magnitude, as 
forecast by (2. 13). 



3 Orders of magnitude 

The transformation process in section 2 is carried out for m = 1, 2, M — 1 and, as in j^EI, it can be 
formalised as an algorithm if the various o— estimates are expressed more precisely in terms of suitable 
orders of magnitude. Here we consider the situation where the orders of magnitude involve powers of x~ a 
with a given constant a(> 0). More general schemes are also possible, as in ffl, section 3 and example 
5.2], but we shall not elaborate further here. 



We suppose, then, that ^ = M — 1 in ( 2. 5) and 

V-i = o(x~n (1 < j < M - 1) (3. 1) 

in (2. 6), while the prescribed accuracy for E\ in ( p. 7| ) is stated as 



Ex{x) =0(x- Ma ). (3.2) 
In practice, ( |3. l[ ) will be an exact order of magnitude except that we are now allowing the possibility 



that some Vji with that exact order may be absent from (2. 5). We make similar order assumptions 
concerning p, A and Ai in ( |2. l[ ) and ( [2. 2| ): 

1/ P (x) = 0{x l ~ Ka ), (3. 3) 

\/\{x) = 0(x- La ), (3. 4) 

Ai(ic) = Di +0{x- a ), (3. 5) 



Ai(a:) =0(x- a ), 



(3. 6) 



where K and L are positive integers and D\ is constant. There is a further set of assumptions concerning 
derivatives to state, necessitated by the presence of P' m in (2. 14) and the definition of P rn in ( |2. 15| ). 
We assume that the O— estimates ( ft. l| ) and ( |3. 3| ) - ( |3. 6| ) can be differentiated the relevant number of 
times in the obvious way: that is, each differentiation reduces by 1 the power of x on the right - hand 
side. The relevant number in this context is 



M-j for Vji, M-l for A -1 , 
M-2 for p-^A^Ax. 



A simple induction argument, as in 0, section 5], then gives 



V jm (x) = 0(x-( m+ >-V a ) 



(3. 7) 



in ( |2. lip , where now 

fj, = M - m (3. 8) 

and E m = 0(x~ Ma ), being the prescribed accuracy ( 3. 2 ). In addition, by ( 2. 15 ), ( 3. 3 ) and ( |3. 7| ), we 
have 

p- x P' m = 0(x- (m+K ^), (3. 9) 



so that ( [2. 22 ) is satisfied. Further, since A m and A m also satisfy ( |3. 5 ) and ( 3. 6 ), other estimates 
which follow when (3. 7) is substituted into ( [2. 17| ) - (2. 21) are 



U m = 0(x-(™ +i ) a ), f m = 0(.x-( m+i ) a ), 



f m = 0(x- (m+i+1)a ), T m = 0{x^ m+1 ^ a ) 



We note that, when L = 1, U m ,T m and T m all have the same order of magnitude and can therefore be 
combined into a single term U m say. Then the right-hand side of ( |2. 16 ) can be expressed simply as 



(3. 10) 

where U m — 0(x~( m+1 * >a ), and T m — 0(x~( m+2 ) a ). The order estimates now established are the basis of 
the algorithm which formalises the transformation process of section 2. 

Before moving on to the formulation of the algorithm, we end this section by stating the Levinson 
asymptotic theorem as applied to the final system 

Z' M = p(A M + Rm)Z m (3. 11) 

in the process (|2. 9| )-( 2. 10 ). With <fr M as in ( 2. 12 ), we write 



$ Af =dg(<5i, ...,£„). 

Then, by ( |2. 1%) and (|2. 2|) , the theorem 0, Theorem 1.3.1] states that ( |3. ll| ) has solutions 



Zmk (1 < k < n) such that, as x — > oo, 

Z Mk {x) = {e k + m (x)} exp( I p(t){X(t)d k + 5 k (t)}dt) (1 < k < N) (3.12) 



x 

and X(t)d k is replaced by d k in the integral when N + 1 < k < n. Here e k is the unit coordinate vector 
in the fc-direction and rj k {x) — > as x — > oo. The standard conditions in the theorem are that pRm is 
L(X, oo) together with the dichotomy condition, which we can take in the form j|, (1.3.13)]. 

Re .F(a;)has constant sign (either < Oor > 0) (3. 13) 

in [X, oo), where F(x) is any one of 

pMdj-dJ + Sj-Sk] (l<j,k<N) 



p(dj -d k + Sj - S k ) (N+l< j, k < n) 
p{\dj - d k + Sj -S k ) (l<j<N,N + l<k< n). 

If, in addition, 

/>oo 

/ F(x)dx |= oo (j ^ k) 

then, as in |^|, Lemma 3.1], we have 

n QQ P OO 

\Vk(x)\<( \ P {t)\\\R M (t)\\dt)/{l-n \p(t) HI R M (t) || dt) 



x 



X 



(3. 14) 



(3. 15) 



in [X, oo), where || Rm \\— max | \ (1 < i,j < n). Thus, transforming back to the original system 
( p. 1| ) by means of ( p. 9 ), we see that ( |2. 1| ) has solutions of the form (3. 12) but pre-multiplied by 



and the accuracy achieved in Rm is reflected in r\k by (3. 15). We also note that an obvious situation 



where (3. 13) and (3. 14) hold is given by real p, A and dj in conjunction with (3. 3) and (3. 4) 



4 The algorithm for large A 

We are now in a position to use the analysis in the previous sections to develop an algorithm and computer 
code to evaluate the solutions of < \l. l| ). The general strategy that we adopt is broadly similar to that in 
; i.e. we use the asymptotic analysis in the previous sections to develop an algorithm and computer code 
that will evaluate the solution set in the interval [X, oo), < X, and then use a numerical ODE solver 
to compute the solution over [0, X] with initial conditions at X provided by the asymptotic algorithm. 
As in |^], a feature of the method is that X is small, about 10 or 20. However the absolute error in the 
solution at X is also small: in the examples that we report on it is less than 10~ 6 . 
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We first discuss the algorithm to compute the asymptotic solution of t \l. We assume that M — 1 
iterations of the algorithm are to be performed. As in Q the procedure consists of 3 distinct parts which 
we call Algorithms 4.1-4.3. Algorithm 4.1 is devoted to obtaining a set of recurrence formulae which 
define Sj, j = 1, M, in terms of quantities with a lower subscript, the only assumption being that the 
quantities concerned satisfy non-commutative multiplication. However, the extra structure that ( |T 3[ ) 
imposes must be taken into account in developing the algorithm. Also we can no longer assume that 
Vji (1 < j < /i) or E\ are initially zero and this extra complexity must be taken into account. Further, 
the approach to estimating the absolute error has been modified. Instead of, as in Algorithm 6.1 of 
[|| , computing a cumulative total error Ej , which is the total error committed after j iterations of the 
procedure, we denote by Ej the error computed as a consequence of the j iteration. The total error is 
computed from ( 4. 2| ). 



Algorithm 4.1 ( a ) Input M to specify the number of iterations and hence the accuracy (3. <). 



( b ) Initialise D\, Vji (1 < j < /i) and E\ with initial values that define the problem and also define the 
magnitude of the error. 

( c ) For each 

U e {-p^Q'm, -p-'Q'm, U m , T m , f m , f m , V jm (2 < j < /i), V jm Q m , V jm Q m (1 < j < (J,)}, (4. 1) 

determine the least positive integer v with 

P^ 1 U = 0(E 1 ). 

( d ) For r — to v, determine the order cr m +k — fa m + (order ofU) of P^JJ . 
( e ) Update V ktTn+1 = V k , m+1 + {-IfP^U, E m = E m + A m P^ +1 U. 
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( f ) Output S m +i = Vi. m +i- 

As described in section 2, A m+ i and Vi lTn +i are then obtained from S m +i by moving the diagonal entries 
of SVn+i to A m ; the modified matrix S m +i is then renamed Vi. m +i- The extra structure that the problem 
imposes demands that more quantities appear in the algorithm than in the corresponding algorithm of 

The second algorithm performs the same function as Algorithm 6.2 of 

Algorithm 4.2 The quantities Sj and Ej defined by Algorithm J^.l are realised as n x n matrices. 

The final algorithm is concerned with the evaluation of the sup norm of the total error E. This algorithm 
has some important differences from its counterpart in 

Algorithm 4.3 The total error E is given by 

M 

E = Y J A 3 E 3 (I + V 3 ) (4.2) 

i=i 

where we have written 

I + P 3 = (I + Pi)(I + P 2 )...(I + Pj)- 
In general both Aj and (/ + Pj) are small perturbations of /. However a naive use of the Cauchy-Schwartz 
inequality when evaluating the sup norm of E will induce a factor of n 2 - 7 into each term involved in the 
estimation of the norm. This cost may be considerably reduced by noting that 

A j = (I + V j )- 1 

and using this to obtain a bound for Aj as 

Mill<l+ || VjW /(l — n\\ Pj ||). 
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5 Scalar factor p(x) = x 1 



By J3. SQ , the validity of ( |2. 22j ) depends on having if > 1 in ( |3. 3| ). When X = 0-that is, when 
X/p{x) — 0(x)- those terms in ( [2. 15| ) which involve A continue to satisfy ( 2. 22] ), but those which do 



not involve A contribute only 0(xvL) to the left-hand side of (2. 22), and such terms are no smaller than 
O(vij) . It follows that, in order to eliminate terms of the size of V\ m from (2. 14), the definition (^] 
|l6| ) must be modified to incorporate also p~ 1 P! m as far as entries in the lower right-hand block of P m are 
concerned. 

Let us suppose then, for simplicity, first that p(x) = x^ 1 and second that a typical entry in V\ m 
is CijX~ ma , where Cij is a constant, and i and j lie in the range [N + 1, re]. Then the entry in 

-p~ x P' m + V lm + A P m - P m A 



is zero if, in place of the last of (2. 15), we define 



Pij = Vij/ (dj - di - ma). 



(5. 1) 



The remainder of (|2. 15j ), involving A, is unchanged as is U m in fl2. 16|) and ( |2. 17|) . The definition (|] 
|l]) is valid provided that 

ma ^ dj — di 

for 1 < m < M — 1 and all i and j in [JV + 1, n]. This condition is certainly satisfied if, for example, 



dj — di < 



(5. 2) 



In terms of ( 2. 18| )-( p. 2lQ , the modification ( |5. 1| ) leaves unchanged Q m ,T m and T m , while Q m has 
ma subtracted from the denominators of its entries. The discussion concerning orders of magnitude in 
section 3 continues to apply. 
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In order to obtain a computational algorithm that reflects the analysis of this case we make simple 
modifications to Algorithms 4.1 and 4.2 of the previous section. The components of ( [1. 1| ) which involve 
Q rn are omitted, and in Algorithm 4.2, the definition of P m is modified in order to take ( |5. 1| ) into 
account. Algorithm 4.3 is unchanged. 



6 The generalised hypergeometric equation : an example 

The generalised hypergeometric equation 

m 

y^(x)-x a ^2a j x j y ( - j \x) = (6. 1) 

3=0 

provides a situation where the numerical results obtained from our algorithm can be confirmed indepen- 
dently by the known analytic properties of the solutions of ( |3. 1[ ) . Here a and the cij are constants with 
a > —n and m < n — 1. On the one hand, ( |6~1 ) can be written in the form ( |2. l| ) with p(x) — x 1 
The algorithm in section 5 is applicable, and solutions with the asymptotic form given by the Levinson 
Theorem can be computed back to x = to a prescribed accuracy. On the other hand, the analytic 
theory developed in Q gives exact relationships between solutions with known behaviours for small and 
large x, thereby providing an independent check on the efficiency of our algorithm. 
Here we discuss this matter in terms of the example 

y (3) 0) - x 2 y {2) (x) - xy {1) {x) + y(x) = (6. 2) 



which is the case n = 3, m = 2, a = of (6. 1). We quote from [|| for the detailed results which form the 



basis of our discussion. We first note that 

y (x) = x (6. 3) 
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is an exact solution of (6. 2) and that two other solutions, analytic for all x, have the form 

oc oo 

o o 

where cq = do = 1. 

To express ( 3. 2 ) in the form ( |2. l| ), we define the column vector Y = (y,y^ 1 \y^) T and write ( 6. 2 ) 
in the usual way as Y' — AY, where 



A 



1 
1 

— 1 x x 2 

Then, as explained in we make a series of transformations which, put together, are 

1 
' ! = x- 1 

a; 



111 

x 3 1 -1 
x 3 - 1 2x~ 3 



1 










3a;- 3 


1 





z 








1 





l + 3a;- 3 1 1 
x 2 + 3a;~ 4 a; -1 — x^ 1 
x 4 - x 2a;- 2 



This gives the Z— system 



with 



Z' = x~ l {K + R)Z 



(6. 5) 



(6. 6) 



A = dg(X - 3, 1, -1) 



(6. 7) 
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R = 3 



-(X-l)- 1 2(X 2 -l)- 1 
A{X - I)- 1 +AX- 1 -{X - I)- 1 - GX-^X 2 - l)- 1 

o (x + i)- 1 



and where X — x 3 . In order to obtain the form (2. 5), with ix = 2 for example, we write 



R = ZX- l Cx + 3X- 2 C 2 + E x 



where 





-1 







-1 


2 


Cx = 


8 


-1 


, c 2 = 


4 


-1 







1 







-1 



Ex = 3 



-X^iX-l)- 1 2X- 2 (X 2 -1)- 1 
AX^iX-l)- 1 -X- 2 (X - l)- 1 - 6X- X (X 2 - l)- 1 

o o x^ix + iy 1 

On taking the diagonal terms in Cx over to A, we finally obtain ( |2. l| ) with = x~ x and 

Ai = dg(X - 3 - 3X~\ 1, -1 + 3A- 1 ). 



(6. 8) 



(6. 9) 



In ( p~~2| ) we have X = X = x 3 ,D = (1), D = dg(l, -1), A x = (-3 - 3X" 1 ), Ai = dg(0, 3AT- 1 ). Thus (|] 
4[ ) - ( 3. 6| ) are satisfied with a = 3 and L = 1 while, by ( p. 8 ), ( p. 1| ) is also satisfied. The condition (^] 
2[) is also clearly satisfied. The algorithm in section 5 is therefore applicable. 



The asymptotic solution of (6. 6) is given by the Levinson Theorem and, when the first component 
of ( |6. 5| ) is taken, we obtain solutions yoox and y^ of (6. 2) such that 

.1 



D i ~ x exp(-a; ), y^ ~ x~ 



(6. 10) 
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as x — > oo. These solutions arise from the entries X — 3 and —1 in A. The entry in 1 in A also gives rise 



to a solution asymptotic to x, but this solution is associated with yo(x) in (6. 3). 



6.1 Computational results at x = 



In order to test our algorithm, we shall use it to compute the solution at x = 0. We do this since we 
have an independent check on this value provided by a series expansion. The analysis needed for this is 
derived in the next section. In this subsection we show how our algorithms perform for this example. 
We shall use the theory in the above sections to obtain the solution of ( |6. 2| ) which is asymptotic to 



x . The initial conditions required by part (b) of Algorithm 4.1 are obtained from (3. 8) and are 



H = 2, V u = 3X-\Ci - dgCi), V 2 i = 3X- 2 C 2 , El 



Ai is determined by ( |6. 9 ). Algorithm 4.1 gives 



Si - Vu 

q Q'i 

J2 — 



V 11 Qi+T 1 +f 1 + U 1 + V 21 



while Si is obtained directly from (6. 8). Algorithm 4.2 realises Si and S2 as the 3x3 matrices 



Si 



=4 1 



-3 



£ j 



( 6 \ 

-3 n 6 



S-2 



=M 



36(2+7x 3 ) _ 
x 9 "5 

=§■ 



(6. 11) 
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The formulae ( |2. 15| ) together with the modification ( |5. 1| ) yield explicit expressions for P\ and P2 
respectively. 

The Levinson asymptotic theorem ( 3. 12| ) is applied to the equation ( |2. 10| ) with m — M — 3 and 



the result evaluated at x = 10. By ( |3. 2| ) therefore, our accuracy at 00 is 0(x 9 ). Since we are interested 
in the solution asymptotic to a; -1 , ( 3. 12 ) ( but with k = n = 3), ( |6. 9| ) and ( |6. ll|) give 



Z 33 (x) = C{e 3 + V3 (x)}exp ^jf ^(-l + 3t~ 3 - 3f 



where 773(2;) = 0(x 9 ) by ( 3. 15|) . Here, the arbitrary constant C is introduced and chosen so that Z 3 



^33 

is asymptotic to x~ l to within 0(x~ 9 ). In short 

^33(2:) = {e3 + mix)}^ 1 ex P (~x~ 3 + 
and 

Z 33 (10) = (0,0,0.0999000999) T (6. 12) 



while (4. 2) in this case yields an error estimate 

£(10) = 2.09830422 x 10" 



Pre-multiplying ( 6. 12| ) by (/ + Pi) (I + P2) together with (6. 5) gives a value 



Y(10) = (0.0999600993, -0.009984070, 0.0019920140) T . (6. 13) 

This is then used as the initial conditions to the NAG |7j library numerical ODE solver D02NMF, which 
is used to evaluate 

Y(0) = (1.87778537, -1.76303921, 1.99999920) T (6. 14) 



which is accurate to an error of 10 7 when compared with the result (7. 3) from the next section. 
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We remark that in order to further verify the accuracy of our method we have implemented the 
numerical procedure described above in the interval arithmetic code CXSC. This together with the 
Taylor based ODE solver || has been used to compute enclosures for Y(0) with initial values provided 
by the asymptotic analysis. These enclosures are 

Y(0) = (1.8777?°, -1.763< ffi) (6. 15) 

and contain both the computed solution ( |6. 14| ) and the value ([7. 3|) given by the series in the next 
section. We further remark that although the error in the initial condition ( |6. 13| ) at x — 10 is of order 
1CP 8 , both integrators have lost an order of magnitude in accuracy in integrating over (10,0). 



7 Analytic theory 



We turn now to the analytic theory of ( 6. 2j ) as given in || pp 78-97]. In the notation of || p. 78], we 
have n — 3, p = 2, fix = — /3g = 1, K = | and 8 = — 1 in the case of ( |6. 2| ). We consider first the solution 
^3,2(2;) given by J|, (3.6.4)]. In order to relate this solution to the real- valued solutions ( 5. 3 ) and ( |S. 4| ), 
we define 



V(x) = ±{V 3 , 2 (x) + V 3 Axe 2 ^ 3 )} 



2- 



^ k (3 2 / 3 x) fc cos(7rfc/3) 



fc=0 



Because of the T— functions in the denominator, V(x) is a power series in x 3 with constant leading term 



\2> 3 / 2 , together with a single term in x. Hence, by (6. 3) and (6. 4), 



V{x) = \2? /2 yx{x) - ^^yo(x). 
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The asymptotic form of V(x) is determined by J8|, (3.6.5) - (3.6.6)], and it is sufficient for our purposes 
to note that 

V^)— ^3 3 / 2 x- 3 exp(i* 3 ). 



On comparing with (S. 10) and noting that y^ and j/q are sub-dominant at oo , we conclude that 



2/i (x) = -Vooi{x) + c 1 y oo2 (x) + dxyoix), 



(7. 1) 



where c\ and d\ are constants which may or may not be zero. 

Next we consider the solution Wz,i (1, z) given by || (3.7.4)]. Again, to keep to real-valued solutions, 
we define 

*r(4 + .^ 



W(x) = lw 3 , 2 (l,xe^) = £{-1)* ^yjL (?'*z) k 

K.L{ 3 3 ) 



fc=0 



As in the case of V(x), the T— function in the denominator and the definitions (3. 3) - (6. 4), show that 

(7. 2) 



o4/3 o 

W{x) = Zy x {x) + ^iyJ/ 2 (x) - 3 2 / 3 T(-)y (x). 



The asymptotic form of W(x) is determined by |8|, (3.7.5)] and again it is sufficient for our purposes to 
note that 



W(x) 



ft therefore follows from ( ]6. 3[ ) and ( |6. 1C| ) that 



Then ( 7. 2 ) gives 



3 i/3 



yoo2 (x) = 2.3- 1 / 3 T(^)y 1 (x) + y 2 {x) - 2.S- 2 ^{T(^} 2 y (x). 



In particular, we have 



2/oo2(0) = 2.3~ 1/3 r(-) = 1.87778588, 



(7. 3) 
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in agreement with the first component of (6. 14) and (6. 15) which is the value of y<x>2(0) produced by 
our algorithm. 
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